Vacancy Induced Splitting of Dirac Nodal Point in Graphene 
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We investigate the vacancy effects on quasiparticle band structure of graphene near the Dirac 
point. It is found that each Dirac nodal point splits into two new nodal points due to the coherent 
scattering among vacancies. The splitting energy between the two nodal points is proportional to 
the square root of vacancy concentration. In addition, an extra dispersionless impurity band of zero 
energy due to particle-hole symmetry is found. Our theory offers an excellent explanation to the 
recent experiments. 
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One of the important features of pristine graphene is 
the Dirac nodal structures centered at two inequivalent 
corners K and K' of the first Brillouin zone with lin- 
ear energy dispersion [IJ [2] . The common belief is that 
the Dirac nodal structure is robust against the short- 
ranged potential scattering because the Fermi wavelength 
diverges at the Dirac nodal point and the electronic state 
correction due to the scatters is negligible [3J. However, 
the recent progresses in the classical wave physics demon- 
strates that local resonant structures can dramatically 
modify waves whose wavelengths are several orders of 
magnitude larger than the structure sizes [1H5]. Vacan- 
cies as well as various chemical adsorbates in graphene 
can create resonant states in the vicinity of the Dirac 
point. An effect analogous to the classical wave is ex- 
pected. Those locally resonant states should dramati- 
cally change the graphene electronic structures and trans- 
port properties near the Dirac nodal point. Indeed, the 
angle-resolved photoemission spectroscopy (ARPES) in- 
dicates the opening of a tunable band gap near the Dirac 
point and the formation of a dispersionless impurity band 
in hydrogenated quasi- free-standing graphene 0[5]- An- 
other study of hydrogenated graphene on SiC showed sig- 
nals of a metal-to-insulator transition (MIT) due to the 
electron localization [9j . Away from the charge neutrality 
point, the transport measurements demonstrated a sub- 
linear carrier dependence of the conductivity [TU] . Within 
the Boltzmann transport framework, a sublincar conduc- 
tivity in charge density was predicted |llj . Away from 
the Dirac nodal point, this theoretical prediction agrees 
well with experiments. However, it fails to explain the 
transport behavior near the nodal point, which is sup- 
posed due to the breakdown of Boltzmann transport the- 
ory there [TH[T3]. Nevertheless, those predictions have 
not included the influence from the substantial change 
of electronic structure arising from the resonant scatter- 
ing. In fact, numerical simulations show that a dramatic 
change in the density of states (DOS) occurs near the 
nodal point P3MT3]. Therefore, a deep understanding of 



the resonant scattering effects is needed. 

In this letter, we present a calculation on the effects of 
the vacancy resonant scattering on Dirac nodal structure 
in graphene. The quasiparticle dispersion is extracted 
from the spectral function A(k, E) which can be calcu- 
lated by extending the well developed Lanczos approach 
[Tfi] . In contrast to the previous theoretical studies of the 
spectral function by the average T-matrix approximation 
(ATA) [iTlflTj] or the self-consistent T-matrix approxima- 
tion (SCTA) [T71 [2U] , our proposed method is more gen- 
eral and nonperturbative, including all the coherent mul- 
tiple scattering contributions. The ATA considers only 
the averaged contribution of impurity potentials with- 
out interference effect from different impurities while the 
SCTA includes the partial contribution from quantum in- 
terference. We found that, instead of only one peak in the 
spectral function for a given momentum in the case of a 
weak short-range scattering [TS] , the resonant scattering 
yields multi-peaks in the spectral function. This leads to 
a complete change of quasiparticle band structure. Each 
Dirac nodal point splits into two new nodal points and 
the splitting energy between the two nodal points is pro- 
portional to the square root of vacancy concentration. 
In addition, an extra dispersionless impurity band is de- 
veloped at zero energy. The fact that we can obtain this 
existing band proves the accuracy of our method. Our re- 
sults also suggest that the Boltzmann theory could work 
well near the Dirac point if the modified energy disper- 
sion is taken into account. 

7r-electrons of pristine graphene can be modeled by a 
tight-binding Hamiltonian on a honeycomb lattice of two 
sites per unit cell, Hq = t ^ \i >< j\ + h.c, where t 

<ij> 

is the hopping energy between the nearest neighboring 
atoms. We consider vacancy disorder potential V in this 
work because it gives rise to the resonant scattering which 
may also arise from hydrogen or fluorine adsorbates [3T] . 
Vacancies are introduced by randomly removing lattice 
sites with probability ni mp (vacancy concentration). In 
the presence of a disorder potential V, the single electron 
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properties can be obtained from the ensemble-averaged 
Green function G(k±,E) = < k± | E+i } Ha _ v |k± >, 
where |k± > is eigenstate of Hq. This Green function can 
be obtained numerically by using the Lanczos recursive 
method [T5J [Jj>J I2"2"H!M] . To numerically obtain an exact 
ensemble- averaged Green's function near the Dirac point, 
a large lattice containing millions sites (4800 x 4800) is 
used. The large samples guarantee that the calculated 
Green's function is free from the finite size errors. In the 
following calculations, the broadening parameter is set to 
be r\ = 0.001* JTB]. 

Self-energy function E is defined in the Dyson's equa- 
tion as G(k,E) = G (k, E) + G (k, J5)E(k, E)G(k, E). 
Thus one has S(k,J5) = G„ 1 Ql, E) - G'^k, E) [23126]. 
Fig. 1 is the energy-dependence of calculated real (a) 
and imaginary (b) parts of the self-energy function in 
the conduction band for various momenta k = 0, 0.021 
and 0.042 along K — M direction in the Brillouin zone. 
To have a visible effect of the resonant scattering, a rel- 
atively large vacancy concentration rii mp = 0.1% is used 
in the calculation. Within energy range of [—0.02*, 0.02*] 
(depending on the vacancy concentration) , both i?eE and 
TroE depend explicitly on the momentum, as shown in 
Fig. l(c-d). This momentum dependence of the self- 
energy results in the failure of the ATA or SCTA [27]. It 
also shows that the effective homogeneous medium ap- 
proximation is invalid around the Dirac point. Beyond 
this energy range, however, our simulations show insensi- 
tiveness of the self-energy function to wave vector k, and 
E varies only with the energy. More interestingly, a spike 
in i?e£ around the Dirac point is observed in Fig. 1(c), 
in contrast with the ATA and SCTA [TT] [20]- Another 
notable feature is that the ATA predicts a peak of JmE 
at the Dirac point. Instead, a dip at the Dirac point 
(Fig. 1(d)) is observed. Though the self-energy function 
shows a complicated momentum and energy dependence 
around the Dirac point, we find that it still satisfies the 
Kramers-Kronig relation [35]. In order to have a better 
picture of the quasiparticles of the system, we compute 
the spectral function below. 

The single-particle spectral function relates to the 
Green's function through A(k±,E) = — ImG(k±, E)/ir 
[2l>] . Generally speaking, the spectral function Ao(k, E) 
is a (5-function in the absence of disorders, reflecting that 
the wave vector k is a good quantum number and all 
its weight ratio is precisely at energy E — E^±. In the 
presence of disorders, the <5-peak is broadened due to fi- 
nite life time of quasiparticles, resulting from the disor- 
der scattering effect. The linewidth of the peak is given 
by ImYj(E) that measures the elastic relaxation lifetime 
r e , r e = _2imS(E) 1 These general features are indeed 
observed in the weak scattering cases [16]. The spectral 
function is qualitatively different in the strongly resonant 
scattering regime. Taking A(k = 0+,E) in Fig. 2(a) as 
an example, the spectral function is surprisingly split into 
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FIG. 1: (Color online) Real (a) and imaginary (b) parts of 
self-energy function E as a function of energy E for various 
wave vectors: k = (black), 0.021 (red) and 0.042 (blue) (in 
unit of 1/a). The wave vector varies from the K point to M 
point in the Brillouin zone (Inset of (a)), (c) and (d) are the 
enlarged views of (a) and (b) near the Dirac point {E = 0). 

three peaks: One broadened peak (p_) centering in the 
hole regime (E < 0), one broadened peak (p+) centering 
in the electron regime (E > 0) and one sharp peak (po) 
at charge neutrality point (E = 0). The p_ peak moves 
toward E = as the wave vector k increases while its 
height is reduced and its width increases, po peak posi- 
tion do not change with k while its height decreases with 
k. Meanwhile, the p + peak moves away from E = and 
its height is significantly increased and its width is nar- 
rowed. When the wave vector exceeds a threshold value, 
the p- and the po peaks disappear. 

By tracing the trajectories of the peaks of the density 
plot of spectral function A(k, E) = A(k+ , E) + A(k- , E) 
in the k — E plane, shown in Fig. 2(b), one can obtain 
the dispersion relation E(k) shown in Fig. 2(c). Many 
other physical quantities like the group velocity and elas- 
tic scattering time can be directly extracted from Fig. 
2(b-c). Away from the Dirac point, the p + peak domi- 
nates the spectral function and the effective energy dis- 
persion approaches the linear behavior (red square dotted 
line). Around the Dirac point, Fig. 2(b) clearly shows the 
anomalous band structure of quasiparticle. The disper- 
sionless impurity band (green diamond dotted line) corre- 
sponding to the po peak in the spectral function indicates 
the localized states relating to the kink in i?e£ and the 
dip in JmE. The existence of these localized states are 
originated and protected by the particle-hole symmetry 
28 : Each vacancy deducts one level from the continue 
band, and at the same time pulls out another level from 
the continue band to a highly localized zero energy state 
(hopping is not allowed in order to be dispersionless) 
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FIG. 2: (Color online) (a) Quasiparticle spectral function 
A(k+, E) plotted as a function of energy E for various wave 
vectors at vacancy concentration n imp — 0.1%. A(k—,E) 
can be obtained by reflecting A(k+, E) around E = due 
to the particle-hole symmetry, (b) Density plot of spectral 
function A(k, E) = A(k+, E) + A(k-, E) in k - E plane for 
Tiimp = 0.1%. (c) Dispersion relation (dotted line) extracted 
from (b). The blue circles and green diamonds are the new 
resonant band and dispersionless impurity band, respectively. 
For comparison, the linear dispersion of clean graphene is 
plotted as the black solid lines. 

|29j . Using the exact diagonalization method, one can 
also see from the inverse participation ratio analysis that 
these impurity states are highly localized around each va- 
cancy |14j . Applying the Stoner theory for magnetism, 
this zero- bandwidth dispersionless band may lead to fer- 
romagnetism for the system when the electron-electron 
interaction is switched on. This observation is consis- 
tent with the recent experiment that directly confirms 
the dispersionless impurity band in graphene [8]. More 
importantly, a new dispersion quasiparticle band (blue 
circles) with finite lifetime corresponding to the p_ peak 
in the spectral function is observed in our calculations. 
This indicates that a splitting of host band of pristine 
graphene happens near the Dirac point. 

The physical origin of above results comes from va- 
cancy effect. Firstly, they create zero energy disper- 
sionless impurity band, forming resonance scattering cen- 
ters. The strong resonance scattering around the Dirac 
point where the electronic de Broglie wavelength A is 
much bigger than the average vacancy-to-vacancy dis- 
tance L v ~ 1/ \Jrii mp . The coherent scattering between 
the neighboring vacancies mix K and K' points [30] so 
that a band gap is opened at the Dirac point. This cou- 
pling between the resonant scatters results in the ap- 
pearance of p_ and p+ peaks for both A(k+, E) and 
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FIG. 3: (Color online) (a) Quasiparticle spectral function 
A(k = 0, E) as a function of energy E for several vacancy 
concentrations: n imp = 0.01% (red), n imp = 0.1% (blue), 
and nimp = 0.2% (black). AE is the splitting energy between 
the two resonant peaks, (b) The vacancy concentration de- 
pendence of the splitting energy AE (red dots). The x-axis is 
in (ni mp ) 1//2 . The black line is the linear fit AE = A^imp) 1 ^ 2 
with A = 1.76. 

A(k— , E). In fact, the direct evidence of the correspond- 
ing splitting of spectral function was indeed reported in 
hydrogenated graphene by ARPES [2 [8] . As the wave 
vector exceeds a threshold so that the wavelength is small 
A < L„, the quasiparticle is insensitive to the coupling 
of resonant states between neighboring vacancies in the 
low density. Hence, the new dispersive band disappears 
for the short wavelength. This is why ATA approximate 
approach works well far enough from the Dirac point. 

Coming back to the spectral function at wave vector 
k = (as shown in Fig. 3(a)), one striking feature is the 
splitting energy AE (peak-to-peak distance between p + 
and p- peaks). AE is the splitting energy between two 
nodal points at k = (neglecting the impurity band). 
The vacancy concentration dependence of AE is plotted 
in Fig. 3(b) , where the x-axis is in the square root of rii mp 
that measures the inverse of vacancy-to- vacancy distance. 
The square root dependence of AE oc ^/ni mp oc 1/L V 
supports the picture that the coherent scattering between 
vacancies play an important role. 

The qualitative changes of the original linear disper- 
sion relation around the Dirac point give rise to a com- 
pletely different DOS near the zero energy. The DOS p 
can be directly computed from the spectral function by 
p(E) = ./!(&, 23). As shown in Fig. 4(a), a shape 

singular peak due to the dispersionless impurity band ap- 
pears at E = 0. The smooth part of DOS has also a peak 
near E — that is from the two extra dispersion rela- 
tions (blue circles in Fig. 2c). The DOS far away from 
E = are essentially from the original linear dispersion 
relation of pristine graphene. The width of the energy 
regime influenced by vacancies is approximately propor- 
tional to AE ~ ^n imp . All of these results agree exactly 
with the the previous numerical simulation |14[ I15j and 
recent APRES experiment [S]. Thus, it proves that our 
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FIG. 4: (Color online) (a) Comparison of the DOS obtained 
from our spectral function (red line) with the numerical sim- 
ulation of reference |15| (green dots) for rii mp — 0.1%. Ex- 
cellent agreement demonstrates the accuracy of our calcula- 
tions, (b) Conductivity as a function of charge density n e for 
riimp = 0.01% (blue) and riimp = 0.1% (red), respectively. In- 
set: The conductivity plateau width in terms of energy AE C 
vs splitting energy AE shown in Fig. [3] AE C can be ob- 
tained from the relationship n e — J Ef p(E)dE. Black line is 
the linear fit of AE C = AAE with A = 0.95. 

spectral analysis captures the correct physics near the 
Dirac point. Furthermore, it shows that ATA or SCTA 
fail to describe this interesting and important feature of 
DOS [17]. 

The new feature of DOS near the zero energy has im- 
portant consequences in electron transport. To clearly 
see it, we calculate the conductivity by using the Kubo 
formula at zero temperature, 

° xx = ^f2 Tr \^G R {E) ]x G A {E)\. 

Neglect the corrections to current vertex and use the 
Drude-Boltzmann approximation G R G A ~ G R G A [27] , 
our calculated conductivity as a function of charge den- 
sity for n imp = 0.01% and 0.1% is plotted in Fig. 4(b). 
A plateau of conductivity 4e 2 /nh is observed near the 
charge neutrality point. The width of the plateau in- 
creases with the vacancy concentration. Away from the 
plateau, the conductivity increases with charge density 
sublinearly. These results are consistent with those of 
the numerical real space calculations O EU [35] . Inter- 
estingly, plateau width in terms of energy AE C is almost 
equal to the splitting energy AE as shown in the inset 
of Fig. 4(b), which implies that the origin of the plateau 
is the anomalous band structure near the Dirac point. 
This conductivity plateau provides an interpretation to 
the experimental observation |33) . It should be pointed 
out that our calculations do not include the contribution 
from the impurity-band-induced variable range hopping 
conduction which may also be important in understand- 
ing the experiments at finite temperatures [SJ [TU] . 

In conclusion, we have studied vacancy induced res- 
onant scattering effects on the one-electron properties 
of graphene. From our accurate spectral function, a 
new quasiparticle dispersion band due to the coherent 
scattering between neighbor vacancies is predicted, and 



each Dirac point splits into two new nodal points. Fur- 
thermore, a dispersionless impurity band is developed at 
zero energy. This result sheds light on the characteris- 
tics of electronic structures and transport properties of 
graphene in the presence of resonant impurities. 
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